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Abstract:  A  prototype  expert  system  for  the  treatment  of  stochastic  control  and  nonlinear  signal 
processing  problems  is  described  with  several  illustrative  examples.  The  system  is  written  in 
MACSYMA,  LISP,  and  PROLOG.  It  accepts  user  input  in  natural  language  or  symbolic  form;  it 
carries  out  the  basic  analysis  of  the  user’s  problem  in  symbolic  form  (e.g.,  computing  the  Bellman 
dynamic  programming  equation  for  stochastic  control  problems  or  the  Zakai  equation  and  the 
estimation  Lie  algebra  or  likelihood  ratio  for  nonlinear  filtering  problems);  and  it  produces  output 
in  the  form  of  automatically  generated  FORTRAN  code  for  the  final  numerical  reduction  of  the 
problem.  The  system  also  has  a  module  using  PROLOG  which  can  check  the  well-posedness 
(existence  and  uniqueness)  of  certain  classes  of  linear  and  nonlinear  partial  differential  equations 
specified  in  symbolic  form  by  computing  a  natural  Sobolev  space  for  the  solutions  and  verifying 
classical  existence  and  uniqueness  criteria  for  the  given  equation  using  MACSYMA  for  the  compu¬ 
tations  and  PROLOG  for  the  logical  analysis.  Sample  sessions  with  three  of  the  modules  of  the 
system  are  presented  to  illustrate  its  operation.  The  status  of  the  system  and  plans  for  its  further 
development  are  described. 


1.  Introduction 

In  this  paper  we  shall  describe  a  prototype  expert  system  using  symbolic  manipula¬ 
tion  languages  developed  for  the  analysis  and  design  of  stochastic  control  and  signal  pro¬ 
cessing  systems.  The  system  has  four  major  components:  (i)  a  modular  system  of  pro¬ 
grams  written  in  MACSYMA1  which  “solve”  certain  stochastic  control  and  filtering 
problems  in  symbolic  form;  (ii)  an  “intelligent”  command  language  interface  using 
OBLOGIS,  a  PROLOG  system  written  in  LISP  [1];  (iii)  a  "theorem  proving  module,” 
also  using  OBLOGIS,  capable  of  checking  the  well-posedness  (existence  and  uniqueness 
of  solutions  on  a  Sobolev  space  determined  by  the  module)  of  certain  classes  of  linear 
and  nonlinear  PDE’s  in  symbolic  form;  and  (iv)  a  general  purpose  module  for  generating 
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FORTRAN  programs  from  the  symbolic  manipulation  modules  of  the  system.  The  func¬ 
tion  of  the  system  is  to  accept  input  from  the  user  in  natural  language  with  model  equa¬ 
tions  expressed  in  symbolic  form,  to  “expertly”  (automatically)  select  a  solution  tech¬ 
nique  for  his  control  or  filtering  problem,  to  reduce  his  model  equation  by  symbolic 
manipulations  to  a  form  appropriate  for  the  technique,  checking  well- posed  ness  of  the 
model  along  the  way;  and  to  automatically  generate  a  numerical  language  (FORTRAN) 
program  realizing  the  solution  algorithm. 

In  this  way  the  system  serves  as  a  kind  of  symbolic  compiler  for  control  systems 
engineering.  This  is  its  intended  purpose.  See  Figure  I. 


The  system  which  now  exists  is  capable  of  treating,  among  other  problems2,  sto¬ 
chastic  optimal  control  problems  by  symbolically  evaluating  the  Hamilton-Jacobi 


Figure  1.  Structure  of  the  expert  system. 


2It  can  also  handle  stochastic  approximation,  decentralized  stochastic  control,  optimal  control  by  the 
Maximum  Principle,  system  identification,  and  certain  nonlinear  control  and  filtering  problems.  Several  oth¬ 
er  modules  are  under  development,  including  a  program  for  automatic  implementation  of  the  Hunt-Su-Meyer 
method  for  exact  linearization  of  nonlinear  control  systems  [2]  (suggested  by  C.I.  Byrnes). 


dynamic  programming  equation  for  the  evolution  of  the  optimal  return  function,  select¬ 
ing  from  four  methods  for  the  analysis,  and  then  producing  a  FORTRAN  code  for  the 


numerical  solution  of  this  equation  [3].  To  execute  the  system,  one  inputs3  the  system 
model  equations,  identifying  state  and  control  variables,  and  specifies  in  general  terms 
the  type  of  optimization  problem  involved.  The  cost  function,  the  initial  conditions,  and 
any  other  boundary  conditions  are  also  entered.  The  MACSYMA  portion  of  the  code 
carries  out  the  calculations  involved  in  implementing  dynamic  programming  for  the 
optimization  problem.  That  is,  it  carries  out,  in  symbolic  form,  the  differentiations, 
function  minimizations,  and  routine  algebra  necessary  to  compute  the  Hamilton- Jacobi 
equation  which  describes  the  evolution  of  the  optimal  cost  and  the  functional  form  of  the 
optimal  control  law.  This  portion  of  the  code  involves  a  grammar  based  on  MACSYMA 
expressions.  A  subsequent  portion  of  the  code  which  is  written  in  MACSYMA  and  LISP 
evaluates  FORTRAN  expressions  for  the  appropriate  MACSYMA  expressions  and  pro¬ 
duces  “dimension”,  “write”,  “format”,  “goto”,  “logical  if”,  “numerical  if”,  and  other 
standard  FORTRAN  lines  for  the  construction  of  a  FORTRAN  program.  The  FOR¬ 
TRAN  (subroutine)  is  written  into  a  file  using  a  MACSYMA  command.  For  typical  sto¬ 
chastic  control  problems  execution  of  the  MACSYMA  program  and  generation  of  the 
FORTRAN  code  takes  about  30  seconds.  The  programmer  can  at  this  point  enter 
numerical  values  for  the  system  functions  and  data  and  then  execute  the  FORTRAN 
code.  (These  could  be  also  entered  automatically  from  preassigned  files,  and  the  FOR¬ 
TRAN  code  executed  automatically.)  A  sample  session  with  this  component  of  the  sys¬ 
tem  is  summarized  in  the  next  section. 


The  advantages  of  this  system  are  clear.  Given  a  new  control  problem,  the  time 
involved  in  analyzing  the  dynamic  programming  equation  and  then  writing  the  FOR¬ 
TRAN  code  to  execute  it  is  eliminated.  The  mistakes  and  the  time  required  to  test  and 
debug  the  FORTRAN  code  are  also  eliminated.  This  is  a  major  advantage  of  the  system 
in  its  present  form.  Most  importantly,  the  system  allows  the  engineer  to  interact  with 
the  computer  for  design  in  terms  of  symbolic  expressions.  In  this  way  he  can  modify  his 
analysis  or  design  problem  by  modifying  the  symbolic  functional  form  of  the  model.  The 
FORTRAN  subroutines  that  he  might  have  to  modify  by  hand  to  accomplish  this  in 
conventional  design  procedures  are  written  automatically  for  him.  The  advantages  of 
working  at  the  higher  level  are  clear.  As  we  shall  show  later,  the  system  may  also  be 
used  for  basic  research  in  control  theory. 

The  stochastic  control  portion  of  the  system  also  has  production  rules  written  in 
OBLOGIS  which  take  advantage  of  the  of  the  special  structures  associated  with  such 
problems  to  select  or  guide  the  selection  of  analytical  or  numerical  procedures  to  treat 
specific  applications. 


1.1  Related  work 


Several  systems  for  computer-aided-design  of  control  systems  have  been  developed 
recently: 

(1)  ORACLS  by  NASA  [4]; 


®Using  a  special  grammar  ot  MACSYMA  expressions,  or  a  natural  language  interface.  See  section  2.1. 


(2)  DELIGHT.MIMO  by  the  University  of  California,  Berkeley  and  Imperial  College, 
London  [5); 

(3)  DELIGHT.Maryland  by  the  University  of  Maryland  [6]; 

(4)  CACE-III  by  the  General  Electric  Company  and  RPI  [7]; 

(5)  The  Federated  System  by  the  General  Electric  Company  [8]; 

(6)  and  various  commerical  packages  like  Matrixx,  Program  CC,  and  Cntl-C. 

In  addition,  there  has  been  some  excellent  work  on  general  and  specific  algorithms  for 
numerical  problems  associated  with  control  system  design  [9].  The  cumulative  product 
of  these  efforts  is  a  strong  library  of  numerical  systems  and  packages  for  the  analysis 
and  design  of  numerically  defined  control  systems. 

The  system  in  [7]  is  closest  to  ours  in  spirit.  It  has  an  interface  written  in  LISP 
which  provides  an  “intelligent”  interaction  with  pre-existing  design  packages. 

The  system  described  here  has  a  different  function.  It  is  designed  to  automatically 
generate  numerical  programs  for  the  solution  of  problems  posed  in  symbolic  form.  By 
using  symbolic  manipulation  programs,  the  system  permits  design  efforts  to  take  place  at 
the  “next  level  up”  from  FORTRAN  code.  It  interfaces  more  naturally  with  AI  con¬ 
structs  since  it  is  based  in  LISP  and  involves  symbolic  expressions,  as  opposed  to  the 
numerical  structures  of  FORTRAN.  In  this  way  it  facilitates  the  incorporation  of  an 
intelligent  natural  language  interface.  It  would  be  difficult  to  design  such  an  interface  in 
FORTRAN,  and  it  would  be  equally  awkward  to  write  theorem  proving  modules  in 
FORTRAN. 


1.2  Summary 

In  the  next  section  we  present  a  summary  of  a  sample  session  with  one  module  of 
the  system  with  the  interaction  taking  place  through  the  natural  language  interface.  In 
the  Section  3  we  explain  how  the  MACSYMA  portion  of  the  code  is  developed  for  one 
particular  application  -  the  construction  of  approximate  nonlinear  filters  for  systems 
with  small  parameters.  In  Section  4  we  present  a  sample  session  with  the 
PROOFONPDE  module.  In  Section  5  we  discuss  the  present  status  of  the  overall  sys¬ 
tem  and  plans  for  its  further  development. 

Acknowledgements:  We  would  like  to  thank  C.I.  Byrnes,  P.  Kumar  and  K.  Paul  for 
their  contributions  to  the  University  of  Maryland  portion  of  this  work. 


2.  A  sample  session  for  optimal  stochastic  control 


The  operation  of  the  system  in  treating  a  static,  elliptic  stochastic  optimal  control 
problem  is  illustrated  in  the  transcript  which  follows.  The  problem  is  first  specified 
using  symbolic  expressions.  In  the  interaction,  the  user  specifies  the  states,  controls,  and 
system  structure,  including  the  physical  interpretation  of  the  problem.  In  this  case,  the 
system  is  told  that  the  problem  corresponds  to  the  problem  of  optimal  production  con¬ 
trol  of  a  hydropower  system.  The  states  are  the  volumes  of  water  stored  for  release  in 
two  reservoirs,  and  the  controls  are  the  release  rates.  The  system  retains  this  informa¬ 
tion,  and  returns  it  when  asked  for  the  physical  meaning  of  the  problem.  Notice  that 


the  system  returns  a  summary  of  the  information  it  has  on  the  mathematical  structure 
of  the  problem  at  any  stage,  and  the  information  is  returned  in  the  same  terms  used  to 
introduce  it.  Notice  too  that  the  system  checks  for  well-posedness  of  the  problem  by 
keeping  track  of  dimensions,  indicating  that  the  problem  type  (parabolic  or  static)  has 
not  yet  been  declared,  and  asking  for  the  specification  of  a  discount  factor.  Finally,  and 
most  importantly,  notice  that  the  information  the  user  has  to  enter  corresponds  to  the 
information  necessary  to  write  down  a  dynamic  programming  problem.  The  system  per¬ 
forms  all  the  required  calculations  needed  to  determine  the  Belman  equation  from  the 
system  dynamics  and  derive  the  numerical  implementation  of  the  equation. 

When  the  problem  specification  is  complete,  as  shown  in  the  example,  the  system 
invokes  a  MACSYMA  routine,  in  a  procedure  which  is  invisible  to  the  user,  to  compute 
the  associated  Hamilton-Jacobi-Bellman  (HJB)  dynamic  programming  equation.4  The 
system  will  display  the  HJB  equation  in  symbolic  form  if  asked.  At  this  point,  the  user 
asks  for  the  generation  of  a  FORTRAN  program  to  solve  the  HJB  equation  numerically. 
This  may  be  done  in  a  variety  of  ways.  The  user  may  simply  enter  the  command 
“graph,”  which  is  equivalent  to  asking  the  system  to  compute  the  optimal  return  func¬ 
tion  and  make  a  graph  of  it.  The  system  generates  a  FORTRAN  program  and  tells  the 
user  where  to  find  it;  it  then  proceeds  to  execute  the  FORTRAN  code  and  produce  the 
desired  graph  (shaded  by  hand  in  the  transcript).  The  numerial  results  are  stored  in  a 
file,  and  this  is  also  reported  to  the  user. 

2.1  Interaction  through  the  command  language  interface 

In  the  session  (Figures  2a, b)  the  user  sets  up  the  stochastic  control  problem 
inf  {/0  ’/  (^l(0^2(<).«l(0.«2(t))  dt  |  uvu2} 
f  (u  VU  2,X  ltX  2)  —  U  ?  +  «2  +  X1(l-J1)x2(l-X2) 

+  ((-^-Zi)*2(l-*2))2  +  ((-^ —  a:2)^  1(1  - 

with  the  discount  rate  q  —  1  and  the  state  equations 

dxx(t )  —  [  —  -  u,  ]  dt  +  dw^t) 

dx2{t)  —  {~  -  u2]  dt  +  dw2{t) 

(u  x(t ),«  2(f  ))  e  [0,1]2,  (x  j(<  ),s  2{t ))  e  [0,1]2. 

as  a  model  for  the  control  of  a  pair  of  hydropower  reservoirs  with  states  x1  and  x2,  the 
stored  water  volumes,  and  controls  and  u2,  the  release  rates  of  the  water  from  the 

’in  fact,  the  system  is  capable  of  selecting  from  several  methods  for  treating  a  given  optimization  prob¬ 
lem.  See  section  2.4. 


reservoirs.  The  cost  function  corresponds,  in  effect,  to  the  cost  of  producing  the 
equivalent  thermal  power.  See,  e.g.,  [10].  The  processes  w^t)  and  w2{t)  are  indepen¬ 
dent  standard  Wiener  processes,  and  r ,  is  the  first  exit  time  of  the  state  from  the 
domain  D  —  [0,l]2. 

The  optimality  conditions  for  this  problem  lead  to  the  following  HJB  equation 
-v  +  min{  (uj.Uj)  €  [0,1]2  | 

UX  i  C/X  2 

+  u2  +  ut2  +  (0.5-x1)(l-J2)222  + 

+  (1-^  x)2^  !(0.5-a:2)  +  +  ~ -} 

OX  i  OX  2 

on  the  domain  D  —  {  {xvx^)  €  [0,l]2  }  with  the  boundary  conditions 

v  —  xxx2  for  x1  —  0  and  x2  —  0 

v  =  (l-i  1)(l-x2)  for  Zj  =  1  and  z2  —  1 


In  the  listing  user  input  follows  the  arrow  (==>).  The  system’s  response  is  on 
the  lines  immediately  thereafter.  Note  the  role  of  the  help  facility  and  the  responses  of 
the  system  to  an  incompletely  posed  problem.  The  total  time  required  for  the  session 
was  about  30  seconds  on  the  Honeywell  Multics  system  at  INRIA. 


EYEMPLE  DE  SESS1QK 


hello 

we  shall  try  to  solve  your  stochastic 
control  problem 

do  you  want  to  keep  a  previous  statement? 

no 

do  you  want  to  keep  the  following  session? 

no 

please  give  the  statement  of  the  problem 
(f  you  have  trouble  type  help 

— >Jx1.x2)  Is  a  state  variable  denoting 
the  water-levels 

state  : |x1,  x2) 

physical  meaning  of  1x1 .  x2)  :  water  -  levels 


»'>[ul.u21  is  a  control  variable 
idenotingturbined-water  \ 

control -Hu1.u2l  .  it; 

- . physical  meaning  of  Tut;  u2) :  turblned  - 

-water  -  .  . . ' 

EZ35?^s>g»0i5^5ilhe  drift-^sS^^^®St>- 

B®taihe^irhYrisi6hs»fihe stalYarid  IKeUrirtar- 
BhusTie 

KMiS2“>hsi O.’SaUl .=0 .5hi2  J  is  a  tdrJ Q,®ar 
r*rr  “wra  I  ue  Vfrhs;  I  0 .5 -v .  I  tfjsESfflP  5^2  J  _  1 

7  '--’—value  of  hvio.5.-  u1  .r-'-n--n--0.5..i2'j  -  • 


7  *=>s  *  10.5,0.5]  is Ihe  diffusion 

^diagonal  matrix 

_  2.T..I  diffusion  :  s'  Till!'  7.7 

K;  value  of  i  ?  J 0  5,0 55  ”  :~I  ' 


—  >c 

^»u1  '2+u2‘2+x  t  •  ( 1  -X I  )•  x2* (1  -x2h  =jE- 

r  U0.5-x1)*  x2*  U-x2))‘2+  ’*  J-r 

1  ((0.5-X2)*  x  1  *  ( 1-x  I  ))‘2  is  the  instantaneous 

-  C05t 

instantaneous  :  c 

value  of  c  ;  (0.5  -  xlXl  -  x2px2 

+!  l-x  l)x  H !  -x2)x2r  (I  -  xlF’xUO.S  -  x2) 
4u2^-*U  I  ^ 


•■>what  do  you  know 

instantaneous  .c  ;7'..;''  '*  1 

value  of  c  :(0I5-x1X1  -*2?x2  7  ..  J4 . 

:  +(1  -xDxi 0 -x2)x2_  7i:  :t^±r. 

-4  (1  -  x  1  &  U0J5  .-4(2)4 

±t:tzlvalue  nr'rrifrsrysi 7^ — t"' 

:F7  :  r  drift :  h  :::;V -  , 

tt  ;  ;;  value  of  h  71075  -ju  1  ;0 ,'5:,-JTi2i 
"^T.  '•‘control ;  lul .‘u2l -- 
-physical  meaning  of  Iul,u2l  :tirbined  - 
^■water  •  gr~~  -- 

K  ----  -  -state  :  , 

E2 — physical  -meaning  of  1x1 ,  x2l  :-weter  -4eve's 


■==>the  problem  is  italic? 

-  ihe  problem  is  static  --  V.,as 


»*>dynammicprogramming -equation  2 
a  discount  factor  Is  missing  • "  ”  ''r 

do  you  wish  to  stop  ?  '  ;  _ 


«->q*1  Is  a  discount  factor  -- 
discount  :  q 
value  of  q  : ) 

«*>help 

do  you  need  help  to 

-  give  the  dynamic  of  the  problem  ( 1) 

-  give  the  horizon  and  the  limit  conditions 

(2) 

77  ~-  -  generate  the  bellman  equation  (3) 

*  ‘ . -  generate  the  fortran  program  end  - 

-  numerical  results  (4)  -  --  — *  " 

;  -■  -  ---7-— ssk  for  the  edition  functions  oHhe  data 
MwsVof  the  robot  (5)  -  -  —  .--rV.sJ.-.; 

EiiiSsS-a  please  answerlhe  -correspondihg'fkjmbeh^ 


]^^»ybrmusWri'ciseflfih^^em35maticSb 
Etime7)ependentbylyping:(theprcblem:tslT?r;ir;: 
KtUUcKparaboli  C.Sae3ESESi'5S5SaffiS^BEI 

^Oi»whenihe  problem  jsstatic-you  must  toives- 
®the  discount  factor  .when  the  problem  .*!»«-••••. 
^;7:;iHs  par|b£lic^qutnust5ive  the  horizon  and  / 
gthe  final  cosUjpp-j-  _  , : 

wemr-r  examples :x ,s Us  jt  discount  (f actorJ'^SiP 

Mr-fr^-tt « 1  .Is  the ‘horizon  Jtime.)„»  ’'7'.„77'---  ■ -7 

■i-.  /T~f  Ira  fi nal -cdst^Tt^JTl" 77 T7..777 ' 

7-77'  If  for  instance  the'Jimit  condition  al  7^777'“ 
:77:'  X  «0  is  of  theneuman/dirichlet  type  _*  '77777 
and  equal  to  xr  1  you  type:  "  “•  ; 

x -0  dirichletfneuman  l4x  ;  . 

you  can  substitute  stopping  to  xHrichletTarnd 
reflexion  to  neuman  . - 


Figure  2a.  A  session  with  the  system  for  optimal  stochastic  control 


»ox1  »0  dirichlel  2 

limit  condition  for  x  1  *0  :  1 
dirichlel  type 

“‘delete  limit  condition  *1  ‘0 
limit  condition  deleted 

-**>x1  -0  dirichlel  1  - 

‘“limit  condition  for  xl— 0  -  -  ■ 

•  dirichlet  type  -  -  - 

S-«*‘x2  «0  dirichlet  1 
t  limit  condition  for  >;2 
Edirichlel  type 


■  —  >x  I  —  1  neuman  0  •■■■-= 

Er-limit  condition  for  xl  *1 
wneuman  type ....  - 

^.limit  condition  forj<^«0  : 1- 

gdirlchlet;|yt>fctvj^~f>-.:-.-. .u;_ 

t“‘x2~1  neuman  0  .’••  1 

i  limit  condition  for.>2  *  1  :  0  -= 

neuman  type  ...  .  — . 

.limit  condition  f orixZ^Ojtj - 
pdi r i chi l e tlype^ ~,7ni-ni,Li  JJx  Z 


••‘equation 

hamilton-jacobi  equation  : 

imin((0.5  -  x1)(  i'-.x2)2x2  '  ±  '  f: 

;-  + (1 -x1)xl  (1 -x2)x2 

"-1+ (1-  xlPx  1(0.5  -x2Hut24u22J^1~^. 

Vdv  V  ,  ..d^ 

ij^i0.5'rAl2)---ri6;5 


£v*Horxl->"0 

E-v*H'or-x2-0 


fe-'dv< 


--^-^Oforxi  *1 


BSUx-l  - 

tPZrr-  *D  forx2-l "i~:\ 
»~;dx2  -.  -.V 


:  -do  you  wish  to  stop  ?  .  V 

no  .  ..  .  .  :  ..  s 

-  ■= ‘generate .-m-v v/ 

you  will  find  the  subroutine  in  the  segment 
belman.fortran 

do  you  wish  to  stop  ? 

no 


“‘principal  program 

you  will  find  the  principal  program  in  the 
segment  pp. fortran 

do  you  wish  to  stop  ? 
no 


“‘execute 

«-■  you  will. find  the  fortran  program  in  the 

if  “segment  num. fortran.  .. 

you  will  find  the  numerical  result5-inTile08 
_  ’  where  v  denotes  the  bellman 

'.  function  and  u  the  optimal  control  t?-.-.  ... 
■—  do  you  wish  to  .stop  ?  i;  ■ 

■  no 


Figure  2b.  Continuation  of  the  session 


2.2  The  MACSYMA  grammar  for  the  example 


The  class  of  stochastic  control  problems  treated  by  the  system  is  defined  internally 
by  a  special  grammar  shown  in  Figure  3.  The  grammar  uses  the  Backus-Naur  notation 
where  “|”  stands  for  “or”  and  <word>  stands  for  a  nonterminal  word.  There  are  also 
semantic  constraints  which  define  a  sublanguage  of  this  grammar.  These  constraints 
effect  the  use  of  <y>.  We  shall  not  discuss  this  feature  here. 


<stochastic  -control  -problem  >::=<  domain  >,<inside  -condition  >  ,<boundary -condition  > 

<  domain  >::=  [0.1]'  |  [0,1]"  X  [0,T]. 

<boundary  -conditions  >::=  J2  <boundary  -conditions  >,<  boundary  -element  > 

<Ho*ni*ry -figment  > 

<  boundary  -element  >::=  {<y  >:  <y  >  6  <domain  >,  <x,  >  —  l  }  | 

{<j/>,  <y  >  £  <domain  >,  <x,  >  =  0  | 


{(z.l),  z  €  [0,1]"} 


<x,  >::=  X\  |  z2 


<y  >::=  x  |  (z  ,t ) 

dV 

<boundary -condition  V  =  f  [  —  =  f 

an 

<  operator  >  ::=<  evolution -operator  >  +  <space  -operator  > 

d  V 

<  evolution  -operator  >  ::=  ...  |  — — 

at 

<  operators  >::=  <  operator  >  |  <operator>  +  <operators> 

|  min  {<  operator  >.<  operators  >} 

<inside  -condition  >::==  < operators  >  =  o 

*  qv  <92  V 

<  space  -operator  >:=  -\{<y>)V  +  X]M<  y  >)  - —  +  X)  °‘  ( <  V  > )  "Z- a  +  ('(<y>) 

,=i  9x,  ,  dx/ 

’  BV 

|  min  (tz  £  IRm  :  -  \(<y  >),u)V  +  X)[t,(<y  — 

i-i  ®x< 


+  a,«y  +  C(<y>.«)} 


Figure  3.  Prototype  grammar  for  stochastic  control. 


The  MACSYMA  program  which  sets  up  the  HJB  equation  in  the  example,  written 
in  the  special  grammar  for  stochastic  control  problems,  is  shown  in  Figure  4.  The  gram¬ 
mar  is  highly  efficient,  and  so  the  MACSYMA  program  is  quite  compact. 


2.3  Specialized  natural  language  interface 

The  interface  illustrated  in  the  interactive  session  is  a  program  written  in  LOGIS 
[l l]  which  recognizes  the  specialized  terminology  of  stochastic  control.  The  sentences 
the  user  enters  are  analyzed  and  matched  against  aptterns  of  typical  sentences.  Depend¬ 
ing  on  the  type  of  sentence,  various  actions  are  initiated:  answer  the  user;  access  a  data 
base;  etc. 

For  example,  when  the  user  enters 


let  (  x  v  ]  be  a  state  variable  denoting  water  level 


the  corresponding  internal  pattern  is: 


?-  let  !x  !-  !-  itype  ?-  denoting  !-  Iphysical-meaning 


where  ?-  matches  anything,  !-  matches  one  word,  and  !x  matches  one  word  and  binds  it 


appel():— ( 

fff:xl*(l-xl)*x2*(l-x2), 

fsc:((0.5-xl)*x2*(l-x2))**2+((0.5-x2)*xl*(l-xl))**2-f 
2*x2*(l-x2)+2*xl*(l-xl), 
condi:  [[dirich,  xl*x2],  [dirich,  xl*x2]], 
condim:  [[dirich,  (l-xl)*(l-x2)],  [dirich  ,(l-xl)*(l-x2)]], 
type:  [2,  stat,ecriture,pasmoy,condlim,  condi,  condim,elp,l,difu,  derive, 
plus,dif,[l,l],belm,2,newto,an,ul*pl4- 
u2*p2+ul**2+u2**2+fsc+fff,param,[]], 
hjb(type,prodyn,belman)); 


Figure  4.  MACSYMA  program  for  stochastic  control 


to  x.  The  result  of  the  pattern  matching  is  an  association  list: 


(x  .  \xy,x2],  type  .  state,  physical-meaning  .  water  level) 


The  associated  actions  are: 

•  transform  the  data  and  store  them 

•  verify  the  properties  and  meaning  of  the  data 

•  answer  the  user. 

The  answer  is: 

state:  [x  v  x2] 

physical  meaning  of  \xv  x2\-  water  level 


Through  the  interface  the  user  can  delete  a  variable,  ask  for  the  meaning  of  a  vari¬ 
able,  and  ask  for  all  the  data  and  variables  in  the  system  at  any  point  in  the  interaction. 
The  user  can  request  certain  actions  including  the  generation  of  the  dynamic  program¬ 
ming  equation  associated  with  a  problem,  or  the  generation  of  FORTRAN  code  to  imple¬ 
ment  a  specific  algorithm  for  the  optimization  problem. 

2.4  Rules  for  selection  of  solution  method 

Using  the  logic  programming  facilities  of  the  OBLOGIS  system,  it  is  possible  to 
develop  rules  for  evaluating  the  problem  statements  provided  by  the  user  to  test  whether 
or  not  they  are  complete  and  to  request  additional  information.  In  effect,  the  system  is 
able  to  determine  when  an  optimal  stochastic  control  problem  is  completely  posed.5 
Rules  have  also  been  incorporated  into  the  system  to  select  among  four  different 
methods  for  treating  optimal  stochastic  control  problems.  The  MACSYMA  portion  of 
the  program  is  capable  of  treating  stochastic  optimal  control  problems  using  four 
different  methods:  (i)  direct  numerical  solution  of  the  Bellman  equation  for  systems  with 
small  state  dimension  using  methods  in  [12,13,14];  (ii)  a  decoupling  method  for  systems 
which  have  decoupled  dynamics  and  product  form  objective  function  [15];  (iii)  the  sto¬ 
chastic  gradient  method  [16];  and  (iv)  a  regular  perturbation  method  [17,18]. 

The  problem  statement  supplied  by  the  user  is  tested  to  determine  which  method  is 
applicable.  For  example,  the  (regular)  perturbation  method  is  applicable  if 

•  the  problem  is  a  control  problem 

•  the  time  horizon  is  finite 

•  the  noise  intensity  is  small 

•  there  are  no  state  constraints 

•  there  are  not  control  contraints 


6See  section  4  where  the  issue  of  well-posedness  is  discussed. 


•  the  second  order  differential  operator  in  the  Bellman  equation  is  regular 


These  conditions  are  realized  in  the  following  LOGIS  PROLOG  clauses: 

(+—  (methode  Iperturbation) 

(type  Icommande) 

(horizon  Ifini) 

(petit-bruit) 

(pas-constrainte  letat) 

(pas-constraint  Icommande) 

(non-degenere  lhamiltonien  Icommande)) 

Here  ”  signifies  installation  of  a  clause  and  “!”  introduces  a  PROLOG  constant.  The 
first  list  in  the  clause  is  the  conclusion  and  the  subsequent  lists  are  the  hypotheses  which 
must  be  tested  to  satisfy  the  conclusion.  This  order  reflects  the  backward  chaining 
inference  structure  of  PROLOG. 

2.5  FORTRAN  program  output 

Numerical  evaluation  of  the  solution  of  the  Bellman  equation  is  done  by  a  FOR¬ 
TRAN  program  written  automatically  by  the  system.  A  listing  is  given  in  the  Appendix. 
Note  that  it  is  a  fully  commented,  “stand  alone"  program  written  in  FORTRAN  IV.  It 
contains  some  120  lines,  including  comments.  It  is  produced  in  1-2  seconds  by  the  sys¬ 
tem  in  the  course  of  the  interactive  session. 


3.  Algebraic  and  asymptotic  simplification  of  nonlinear  filters 

The  problem  of  estimating  one  diffusion  Markov  process  {  xt  }  in  terms  of  observa¬ 
tions  of  another  Yt  —  a{  ye ,  s  <<  }  based  on  the  model  (Ito  calculus) 

dxt  —  f  {xt)dt  +  g(xt)dw,  (l) 

dx it  =  h{xt)dt  +  dut ,  0  <  f  <  r  <  oo 


is  fundamental  in  many  engineering  applications.  Here  /  ,  g  ,  and  h  are  smooth  func¬ 
tions,  wt ,  ut  are  independent  standard  Wiener  processes,  and  x0  is  a  random  variable 
independent  of  {wt  ,  ut )  for  all  t .  The  problem  of  recursively  estimating  xt  given  the  <r- 
algebra  Yt  generated  by  the  observations,  may  be  formulated  in  terms  of  the  (Stratono- 
vich)  stochastic  partial  differential  equation  for  the  unnormalized  conditional  density 


du  (t  ,x)  —  (L  *  u)(t  ,x)dt  +  h  {x)u  (t  ,x)dyt 


Lf  u  =  ~  dXI(g\x)  u  )  -  dx{f  (x)  u)  -  ~  h\x)  u  (2) 

u  (0, x  )  =  p0  (x  ),  the  density  of  x0 

That  is,  if  (2)  has  a  nice  solution,  then  the  conditional  density  of  xt  given  Yt  is 

p  (t  ,x)  —  u(t  ,x  )/[ fu(t,z  )dz  ]. 


A  considerable  effort  has  been  devoted  to  the  search  for  (recursive)  finite  dimen¬ 
sional  “representations”  of  estimators  in  terms  of  various  “sufficient  statistics”  of  either 
the  solution  of  (2)  or  other  conditional  statistics  (such  as  the  conditional  mean).  This 
effort,  which  has  produced  few  such  estimators,  has  nevertheless  increased  our  under¬ 
standing  of  the  system  (2)  and  of  the  tools  available  for  its  treatment  [19],  In  [20]  a  sys¬ 
tematic  procedure  was  developed  for  constructing  approximations  to  the  solutions  of  (2) 
in  terms  of  certain  Lie  algebras  of  differential  operators  associated  with  the  equation. 
See  also  [21 ,22,23].  The  procedure  is  based  on  the  presence  of  a  small  parameter  in  the 
model  (1),  and  it  uses  the  methods  of  asymptotic  analysis  to  derive  the  approximations. 

3.1  Asymptotic  analysis  and  the  estimation  algebra 

The  idea  is  as  follows:  Consider  the  simple  perturbation  of  the  Kalman  filtering 
problem  treated  in  [20], 

dxt  =  aXj  dt  +  dwt 

dyte  =  [  xt  +  e{xt)k  }dt  +  dut  (3) 

y0f  =0,  x0  has  Gaussian  density  p0  (x ),  k  >  2. 


The  conditional  density  of  xte  given  Y*  satisfies 

du  €(t  ,x)  =  [L*  -  —  (a:  +  ex  k  )2]u  e(f  ,x  )dt  -f  ( x  +  e  xk  )  u\t  ,x)  dyf  (4) 
2 

u  \0,x  )  =  p0  (x  ) 

L*  u  =  \-dxxu  -  a  dx  {xu  )  . 


Assume  that  u(  has  an  asymptotic  expansion 

u£  =  u0  +  e u  x  +  ezu2  +  ■  •  ■  (5) 

Substituting  this  into  (4)  and  equating  coefficients  of  like  powers  of  e  gives 


(6) 


duQ  —  L*  «0  dt 


+  *  «o  dy,1 


«o  =  Po(x) 


duk  —  L*  uk  dt  +  x  ukdy(t  +  x 1  uk\ 
uk  ( 0,x  )  ==  0,  k  =  1,2,... 


Defining  the  vector 


Un(t  ,x)  = 


«o 

ui 

«» 


we  have 

dUn  (t  ,x)  —  [  (L  *  -  -x^E,  -  xk+1E2  -  — x2kE 3j  U„(t,x)dt 

2  2 

+  [xEl  +  xk  E2]  Un(t  ,x)  dyte 


(7) 


(8) 


where  the  ( n  +1)  x  (n  +l)  matrix  E ,  has  zero  entries  everywhere  except  for  1'  s  on  the 
subdiagonal  (i,l)(i+l,2),...,(n+l,n+2-t)  (i.e.,  on  a  particular  subdiagonal).  Note  that 
Ex  is  the  identity  matrix. 

The  equation  for  u0  is  the  Kalman  conditional  density  equation  (given  t/8e,  s  <  t). 
Therefore,  it  has  a  closed  form  analytical  solution  -  a  Gaussian  density.  This  density 
defines  the  Green’s  function  for  all  the  subsequent  equations,  k  —  1,2,...,  in  (6).  Since 
ttx  involves  this  Green’s  function  and  xku0dytt  as  a  forcing  function,  can  also  be 
computed  explicitly.  It  follows  that  all  the  terms  u1(u2,  •  •  ■  can  be  computed  expli¬ 
citly;  however,  the  complexity  of  the  calculation  increases  very  rapidly.  By  computing 
Un(t  ,x )  we  produce  an  approximation  to  the  conditional  density  which  is  accurate  to 
0(en+1). 

All  these  calculations  can  be  done  automatically  in  symbolic  form  by  MACSYMA. 
To  do  this,  we  use  Lie  theory.  First,  we  rewrite  (8)  abstractly  as 

=  AU„(t)  +  BUn  y  i 

A  —  [{L  *  -  -  xk+1E2  -  l-x^E,}  (9) 

B  —  [x  E j  +  xk  E2] 

Note  that  A  is  a  second  order  differential  operator  and  B  is  multiplication  by  a  matrix 
valued  function.  The  Lie  algebra  generated  by  A  and  B  consists  of  all  linear  combina¬ 
tions  of  A  and  B  and  their  consecutive  commutator  products,  e.g., 


[A  ,B]  —  AB  -  BA 


(10) 


and  subsequent  repeated  products  (taken  as  differential  operators).  This  is  always  a 
finite  dimensional  solvable 6  (e.g.,  “lower  triangular”  structure)  Lie  algebra  for  the  class  of 
filtering  problems  with  polynominal  coefficients  described  above. 

Consequently,  we  can  use  the  Wei-Norman  theory  [24]  to  solve  (9)  globally  in  the 

form 

Un(t,x)  =  [e,l('Mle,2(,M2-  e"(,)A'  ](x )  (11) 

where  (A  ltA  2,...,Ad )  is  a  basis  for  the  Lie  algebra  generated  by  A  and  B ,  and  the  <7,  (0 
are  defined  by  a  lower  triangular  system  of  nonlinear  (stochastic)  ordinary  differential 
equations.  These  equations  must  be  solved  numerically  (in  general).  The  gi(t)'  s  may 
be  thought  of  as  a  finite  dimensional  family  of  (approximate)  conditional  statistics  for 
the  original  nonlinear  filtering  problem. 

3.2  Symbolic  algorithm  structure 

To  implement  the  approximate  nonlinear  filter  using  this  methodology  and 
MACSYMA,  one  must  program  the  following  six  step  procedure: 

(Ml)  Compute  the  differential  operators  A  and  B  in  (9)  in  terms  of  the  original  func¬ 
tions  in  the  filtering  problem  model  (1),  modulo  the  appropriate  power  of  the  small 
parameter  e,  if  one  is  present  in  the  model. 

(M2)  Compute  basis  elements  A  t,A  . . Aj  for  the  Lie  algebra  generated  by  A  and  B  . 

(M3)  Compute  ODE’s  for  the  gt . gd  functions  in  (11)  by  substituting  (11)  into  (9), 

differentiating,  and  inverting  the  resulting  matrix  equation. 

(Fl)  Solve  the  ODE’s  for  the  g> '  s  numerically  using  an  appropriate  numerical  integra¬ 
tion  scheme  for  stochastic  ODE’s  [25]. 

(M4)  Compute  the  approximate  conditional  mean  x  (£  (the  best  mean  square  estimator  of 
xte  given  Ytc)  by  integrating  the  conditional  density. 

(F2)  Represent  the  conditional  mean  computed  in  (M4)  in  FORTRAN  code  in  terms  of 
the  actual  measurements  {  y/,  s  <  t  }  as  input  data  to  the  code. 

The  steps  labeled  (M1)-(M4)  are  done  in  the  MACSYMA  part  of  the  code.  The 
step  (Fl)  is  best  carried  out  in  FORTRAN;  and  so,  the  MACSYMA  step  (M3)  which 
computes  the  symbolic  ODE’s  has  an  auxiliary  component  which  writes  the  FORTRAN 
code  to  numerically  integrate  the  ODE’s. 

The  result  of  step  (F2)  is  the  output  of  the  system  -  that  is,  a  FORTRAN  code  to 
numerically  evaluate  the  best  estimate  of  the  state  in  terms  of  numerical  measurement 
data.  All  the  steps  (Ml)  -  (M4)  and  (Fl)  (F2)  are  carried  out  automatically  after  the 
system  model  (the  functions  /  ,  g  ,  and  h  in  (2))  are  specified. 


6 A  Lie  algebra  is  solvable  if  the  derived  series  of  ideals  L*°*  =  L  ;  Z,("+1)  =  [  L^  \  n  >  0  is  the 

trivial  ideal  {  0  }  for  Borne  n. 


The  system  also  includes  a  test  for  well-posedness  which  consists  of  checking  rela¬ 
tive  growth  conditions  on  the  coefficients  f  g,  and  h  for  existence  and  uniqueness  of  a 
(classical)  solution  to  the  robust  form  of  the  Zakai  equation.  The  criteria  used  are  those 
developed  in  the  paper  [26],  As  is  the  case  in  the  other  components  of  the  system,  com¬ 
plicated  analysis  is  required  to  setup  the  existence  and  uniqueness  test.  This  analysis  is 
done  “automatically”  by  the  program.  Since  the  existence  and  uniqueness  theory  for 
nonlinear  filtering  problems  is  not  complete,  the  program  cannot  give  a  definitive 
analysis  of  a  large  class  of  problems.  Rather,  it  uses  the  classical  existence  and  unique¬ 
ness  theorems  for  parabolic  partial  differential  equations  as  adapted  to  a  class  of  filtering 
models  with  smooth  coefficient  functions.  It  is  not  unlikely  that  more  sophisticated  tech¬ 
niques  could  be  accomodated  by  simple  modifications  of  the  program. 

It  is  possible  to  have  MACSYMA  offer  a  choice  of  numerical  integration  procedures 
to  be  (automatically)  coded  for  computation  of  the  state  estimate.  Our  system  does  not 
now  have  this  capability.  However,  it  does  have  a  second  method  for  computing  the 
conditional  density  based  on  evaluation  of  the  likelihood  ratio  using  formulas  for  integra¬ 
tion  of  stochastic  function  space  integrals  developed  in  [27,28],  The  system  also  has  the 
capability  to  test  the  well  posedness  of  (robust  form  of  the)  Zakai  equation.  In  conduct¬ 
ing  this  test,  using  the  theorems  in  [26],  the  system  tests  the  coefficients  in  the  Zakai 
equation  to  see  if  they  are  polynominals.  If  so,  a  simple  test  suffices  for  well-posedness. 
If  not,  then  a  more  complex  computation  is  executed.  See  [29]  for  details. 


3.3  A  sample  session  with  the  filtering  module 

The  weak  quadratic  sensor  treated  in  the  paper  [20]  is  one  of  the  few  examples  in 
which  the  estimation  algebra  and  the  Wei-Norman  series  representation  of  the  condi¬ 
tional  density  have  been  worked  out  by  hand.  It  is  necessary  to  have  such  examples 
when  programming  in  MACSYMA  to  validate  the  program. 

In  the  session  shown  in  Figure  5  the  program  is  asked  to  analyze  the  simple  filter¬ 
ing  problem 

dxt  —  dwt  (12) 

dyt  =  +  e  xtz  +  dut 

The  system  asks  the  user  to  declare  the  system  as  scalar  or  vector  and  then  to  enter  the 
functions  /  ,g  ,  and  h  in  a  standard  format.  It  then  sets  up  the  Zakai  equation, 
identifies  the  two  operators  in  the  equation  which  will  be  used  to  generate  the  estimation 
algebra,  and  proceeds  to  find  the  representation  of  the  conditional  density  to  the 
appropriate  power  of  e,  in  the  form  of  a  Wei-Norman  series.  The  system  returns  a  basis 
for  the  estimation  algebra  (modulo  the  designated  power  of  e)  and  various  intermediate 
expressions  of  interest,  including  the  matrix  of  structure  constants  for  the  estimation 
algebra.  Its  “final  result”  is  a  set  of  stochastic  ordinary  differential  equations  for  the 
functions  appearing  in  the  exponentials  of  the  Wei-Norman  series. 

The  final  output  of  the  program  is  a  FORTRAN  program  integrating  these  equar 
tions  and  computing  numerical  expressions  for  the  conditional  mean.  The  FORTRAN 


program  written  by  the  system  is  given  in  the  Appendix.  The  complexity  of  the  expres¬ 
sions  in  this  program  effectively  illustrates  the  value  of  a  tool  for  the  automatic  genera¬ 
tion  of  programs.  It  would  be  tedious  to  type  this  code,  and  difficult  to  debug  it  by 
hand.  If  a  single  parameter  in  the  system  model  were  changed,  then  the  entire  code 
might  have  to  be  rewritten.  This  would  be  time  consuming,  and  a  poor  use  of  engineer¬ 
ing  talent.  On  the  other  hand,  the  code  in  the  Appendix  was  generated  in  about  10 
seconds  on  a  VAX  11/785  under  modest  load  conditions. 


4.  A  Sample  Session  with  the  Theorem  Proving  Module 

The  system  contains  a  module  called  PROOFONPDE  for  the  automatic  verification 
of  well-posedness  of  certain  classes  of  nonlinear  PDE’s.  The  module  is  constructed  in 
LOGIS  [11].  It  analyzes  the  coercivity  properties  of  a  given  PDE  and  the  monotonicity 
properties  of  its  nonlinear  components  on  a  Sobolev  space  which  it  constructs  in  the 
course  of  the  analysis.  The  theorems  of  functional  analysis  used  in  the  arguments  are 
encoded  in  rules  in  the  LOGIS  PROLOG  syntax.  The  system  interfaces  completely  with 
MACSYMA;  and  MACSYMA  functions  performing  symbolic  calculations  are  used  when 
necessary.  The  problems  treated  are  in  abstract  form 

Au  =  /  in  Cl  C  IR"  (13) 

BjU  =gj,  j  =0,1  -  1  on  the  boundary  P  of  O 

where  A  and  Bj  are  differential  operators.  The  system  verifies  the  existence  of  a  solu¬ 
tion  to  (13)  and  finds  the  function  space  V  to  which  this  solution  belongs. 

The  function  spaces  used  are  Sobolev  spaces  which  are  represented  internally  as 
MACSYMA  objects.  For  example,  the  space  W™’p  (O)  is  handled  in  the  following  way:  It 
is  typed  by  the  user  in  MACSYMA  syntax  as  space(w,m,p,0)  where  spaceQ  is  a 
MACSYMA  function  that  generates  the  list 

((space  simp)  ((p  m  0  ))  nil) 

which  serves  as  the  internal  representation  of  the  space.  When  the  MACYSMA  disp 
function  operates  on  this  list,  it  prints  W™,v  (n). 

The  main  components  of  the  formal  calculus  (Green’s  function  application,  variar 
tional  formulation,  computation  of  canonical  forms,  etc.)  are  written  in  LISP.  This  is 
done  since  the  main  part  of  the  work  consists  of  manipulating  lists  (MACSYMA  expres¬ 
sions),  and  it  is  most  efficient  to  do  this  in  LISP.  When  necessary,  MACSYMA  functions 
( dijff \  expand,  etc.)  are  called.  The  programming  technique  is  data  driven  programming. 
This  facilitates  the  updating  of  the  knowledge  base. 

The  portion  of  the  system  using  existence  theorems  and  deciding  what  computa¬ 
tions  must  be  done  to  verify  the  hypotheses  is  encoded  in  Horn  clauses  with  the  PRO¬ 
LOG  syntax  of  LOGIS  [11).  This  is  the  core  of  the  system  which  determines  what 
theorem  and  what  method  will  be  applied  in  a  given  case.  I m  example,  the  Lax- 
Millgram  theorem  is  encoded  in  the  following  clause: 


atertOt 


The  model  fo-  the  stochastic  system  in  Jto  for*  is  as  follow* 
dx(t)«  f( s<  t ) >dt  *  g( x( t ) >dw( t ) 
dy  1 1 )«  h<x(t))dt  -  dv(t) 

w(t),v(t)  are  independent,  etenderd  Wiener  processes 
Is  this  the  vector  or  scaler  case? 


input  f(x(tl) 

0; 

Input  g(x(t)) 

t{ 

input  h(x(t>) 

f  •  0 
E  «  1 

2 

h  ■  e  t  ♦  a 

Please  enter  a  number  corresponding  to  the  highest  power  of  'e'  allowed 
1| 

The  Zakal  equation  for  this  problem  is 


c  u 

2  2 

3  u  x  d*  2 

du  «  dt  (»  *  u  i  -  - -  ♦  — )  ♦  dy  (e  u  x  ♦  u  *) 

2  2 


2 

6  u 

2  2 

3  u  x  Ox  2 

pdopa  »  t-  a  v  x  -  — —  ♦  — e  u  x  ♦  u  x) 

2  2 


2 

d  u 

2  2  2 
3  u  *  0*  2  du  du  d  u  du 

kasls  ■  l—  e  u  x  —  -  ■  ■  ♦  — — ,  e  u  x  ,  c  —  x,  u  x,  *— ,  •  ,  e  u  x  e  — 
2  2  dx  dx  2  dx* 

dx 


u,  e  ei J 


t matrix 


l  0  2 

I 

[  1  0 

[ 

l  0  0 

[ 

l  3  0 
t 

l  0  2 

I 

t  0  0 

t 

[  0  0 
t 

t  0  0 

l 

l  0  0 


0  0  0 
0  0  1 
0  1  0 
1  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 


0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  0  0 
0  1  0‘ 
1  0  0 
0  0  0 
0  0  0 


1  3 

3 

o  3 
3 

0  ] 

3 

0  3 
3 

1  3 
3 

0  3 
3 

o  3 
3 

0  3 
3 

o  3 


Exponentiating  the  transition  mat-ix 
Exponentiation  eeapleted 


!»>♦  wymm  or  ».£.•»  kos  tan  HliriliM 


dlddudtlons  >  [0y  cosh(?  t)  »  dy  do.hlt)  -  dy), 

11  dy  dlnhlt)  .  1  dy  »lnh!2  t>),  [dy  co»h(t>),  t-  dy  olnhlt)), 

[dy  co,h(?  t)  -  1  dy  cosh(t)  .  dy),  [-  !  (  dy  dlnh(t)  -  i  dy  co«h!t>), 

1  3 

U  e  dy  ooah  C2  t)  -  ,  dy  »lnh(t)  -  «  ,  dy  oo»h(t)  .  2  ,  dy), 

*  3  t  4 


[-  I  dy  «lnh(t»,  [-  dy  ilnh(2  t)  *  g  dy  coah(2  t)  ♦  2  ,  e  dy  dlnhlt) 
«  *  2  5 


-  g  g  dy  »inh(t) 

3  * 


dy  alnh(t)  ♦  g  g  dy  cosh(t)  -  2 
3  b 


?  2 

g  dy  cosh(t)  -*  g 

*  * 


do  you  want  fortran  generated?  yes;  or  no; 

y«s; 

input  xo 
1} 

input  eps 

U 


To  whet  file  do  you  want  the  fortran  a«nt? 
demfort2; 

(e5)  dsit(); 


Figure  5.  A  sample  session  with  the  nonlinear  filtering  module. 


((solution  space) 

(well  defined) 

(linear  $principal-op  $soiution) 
(elliptic  $principal-op  space) 
(continuous  $principal-op  space)) 


dy)) 


This  is  interpreted  as  follows:  The  problem  (P)  has  a  solution  in  a  space  V  if  (P)  is  well 
defined,  if  the  principal  operator  A  is  linear  with  respect  to  the  solution  u  ,  if  A  is  ellip¬ 
tic  in  V ,  and  if  A  is  continuous  in  V.  Each  clause  triggers  the  execution  of  another 
clause,  MACSYMA  functions,  or  LISP  functions.  All  the  results  or  facts  obtained  during 
the  inference  process  are  kept  in  a  fact  data  base  to  avoid  useless  computations  and 
thereby  to  increase  the  speed  of  the  system. 

In  the  listing  in  Figures  6a, b  inputs  are  on  the  lines  terminated  with  In  the 
first  portion  of  the  program  the  system  explains  its  notation  and  then  asks  the  user  to 
type  in  the  equation  in  a  specific  notation.  It  then  asks  for  the  boundary  condition:  and 
then  it  repeats  the  statement  of  the  problem  in  symbolic  form.  The  system  then  asks 
for  the  nature  of  the  object  “u”,  the  solution.  The  user  supplies  this,  and  the  system 
summarizes  what  it  knows  about  u;  in  particular,  it  does  not  yet  know  the  space  in 
which  u  is  defined  and  it  does  not  know  any  properties  of  u.  This  process  is  repeated  for 
each  element  appearing  in  the  original  equation.  Function  space  and  inequality  proper¬ 
ties  can  be  given  for  each  function. 

When  the  problem  is  completely  specified,  the  system  enters  a  PROLOG  loop.  At 
this  point  the  user  can  ask  several  questions  about  the  problem;  in  the  listing  he  asks  the 
system  to  find  out  whether  or  not  the  problem  posed  has  a  solution.  The  system  prints 
out  everything  it  has  proved:  it  finds  that  the  problem  is  well  defined,  tests  the  linear¬ 
ity,  computes  the  variational  formulation,  and  the  associated  Sobolev  space,  and  verifies 
coercivity.  At  the  end  of  the  manipulations  the  system  reports  that  the  equation  has  a 
solution  in  a  particular  Sobolev  space  which  is  H1( n)  in  the  first  example  and 
H q1  (O)  fl  LT  (n)  in  the  second.  The  total  computation  times  for  the  two  examples  are  9 
and  13  seconds,  respectively. 

The  system  operates  by  putting  the  problem  in  various  canonical  forms  (i.e.,  canon¬ 
ical  forms  of  A  u  and  canonical  forms  for  the  variational  formulation)  where  it  can 
check  the  criteria  spelled  out  in  certain  standard  theorems  for  the  existence  and  unique¬ 
ness  of  solutions  to  nonlinear  PDE’s.  All  the  computations  are  done  in  symbolic  form, 
and  the  conditions  are  stored  and  checked  in  symbolic  form. 


5.  Status  and  objectives  of  the  research 

At  the  time  of  this  writing,  the  work  is  progressing  on  several  fronts.  New  modules 
are  being  added  to  the  system,  including  a  system  for  identification  of  nonlinear 
diffusions  with  jumps,  a  system  for  the  analysis  and  design  of  nonlinear  control  systems 
using  geometric  methods,  together  with  new  optimization  routines,  including  routines  for 
search  and  scheduling  problems.  The  filtering  module  is  being  enhanced  to  handle  vec¬ 
tor  problems  more  efficiently,  and  it  is  being  tested  on  a  variety  of  problems.  The 
natural  language  interface  is  being  redesigned  to  incorporate  more  logic  (PROLOG)  func¬ 
tions.  The  system,  which  now  has  some  15000  lines  of  code  in  LISP,  MACSYMA,  and 
LOGIS,  has  been  ported  on  a  (Symbolics)  Lisp  machine  at  INRIA. 

The  development  of  the  program  so  far  has  concentrated  on  providing  capabilities 
for  executing  specific  control  theoretical  constructions.  In  the  next  phase  of  the  work  it 
will  be  important  to  bring  to  bear  more  of  the  techniques  of  expert  systems  to  provide  a 
systematic  logical  structure  for  the  system.  This  appears  to  be  the  appropriate  long 
term  direction  for  the  work. 


EXEMPLE  1  ! 


(cl)  cogitoO; 

Cogito  :  Version  4  (13/06/84) 
n 

dimension  de  l'espaoe  R  : 

3; 


la  dimension  de  l'espaoe  est  3 

est-ce-correct  7 
oul ; 

equation  a  resoudre  dans  omega  sous  la  forme  "A(u)=f"  : 

-delta(u)+a*u=f j 

conditions  sur  la  frontiere  gamma  sous  la  forme  "(B[ 1 3(u)=0, . . . ,B[m] (u)=0]"  s 
Cdiff (u,nor)=0]j 

le  probleme  a  resoudre  est  ! 

A 

a  u  -  / _ \  u  =  f  dans  omega 

aveo  sur  gamma  : 

du 

dnor 

est-ce-oorrect  7 
ouii 

nature  de  a  : 
ooef f icient ; 

espace  de  a  : 
espaoed.iuf); 

proprietes  de  a  : 
a>1 ; 

fonotion  :  a 


nature  !  coefficient 

espaoe  ;  L(lnf) 

proprietes  :  a  >  1 

est-oe-correot  7 
oui ; 

nature  de  u  : 
solution ; 


fonctlon  :  u 

nature  :  solution 
espace  :  lnconnu 
proprietes  :  nil 
est-ce-correet  7 


fonotion  :  f 


nature  :  sraembre 

espaoe  t  lnconnu 

proprietes  :  nil 

est-ce-correct  7 
oui ; 

=  =  >  ((solution  *espace)) 

**•*  blen_definl  •••» 

****  linealre 
op _princlpal 
solution 

****  formulatlon_variationnelle  •••• 
op_princlpal 
3 

\  [  du  dv  t 

>  I  ( - )  +  I  (a  u  v) 

/  3  dx  dx  ] 

=====  /  1  i  / 

i  =  1 
/ 

[ 

I  (f  v) 

] 

/ 

H(1) 

****  elliptique 

op_prlncipal 

H(1) 

1 

**»»  contlnu  (non  developpe) 

op_princlpal 

HO) 

solution 

HO ) 

==>  fin 

Timex  8809  msec. 


Figure  8a.  Well-posed  ness  of  a  linear  PDE. 


EXEMPLE  2  J 


fonction  i  a 


(cl)  cogitoO; 

Cogito  :  Version  4  (13/06/84) 
n 

dimension  de  l'espace  R  s 
2J 


la  dimension  de  l'espace  est  2 

est-ce-correct  ? 
out ; 

equation  a  resoudre  dans  omega  sous  la  forme  "A(u)=f"  : 
u^abs(u)J‘(p-2)*u-sum(dif  f  (a[i]*dif  f  (u,xCi3)  ,x(i))  ,i,  1  ,n)=f; 

conditions  sur  la  frontiere  gamma  sous  la  forme  "CB[1 3(u)sOt . • . »Btm] (u)«0)*  : 
[usO]; 


le  probleme  a  resoudre  est  : 
n 

\  d  du  p  -  2 

-  >  - (a  — — )  ♦  u  abs(u)  +  u  =  f  dans  omega 

/  dx  i  dx 

*s*s  i  i 

i  =  1 

avec  sur  gamma  : 

U  a  0 

est-ce-correct  ? 
oui ; 

nature  de  u  : 
solution ; 

fonction  :  u 


nature  :  solution 

espace  :  inconnu 

proprietes  :  nil 

est-ce-correct  ? 
oui ; 

nature  de  p  : 
constante; 

proprietes  de  p  : 

p>2; 


fonction  :  p 


nature  :  constante 

espace  :  R 

proprietes  :  p  >  2 

est-ce-correct  ? 
oui ; 

nature  de  a  : 
coefficient; 

espace  de  a  : 

•spaced,  inf ); 

proprietes  de  a  t 
a>i) 


nature  :  coefficient 

espace  :  L(inf) 

proprietes  :  a  >  1 

est-ce-correct  ? 
oui ; 

proprietes  de  f  : 
nil; 

fonction  :  f 


nature  :  smembre 

espace  :  inconnu 

proprietes  :  nil 

est-ce-correct  ? 
oui; 

==>  ((solution  *espace)) 

#,##  bienjjefini  •••• 

***•  bien_defini  (deja  demontre) 

••••  non_lineaire 

op_principal 

solution 


*•••  for»ulation_verletionnelle 
op  principal 

2~ 

*=  =  =  /  /  / 

\  [  du  dv  t  P-2  [ 

>  I  (a - — -)  ♦  1  (u  abs(u)  v)  I  (u  v) 

/  3  i  dx  dx  3  ] 

s===  /  i  i  /  / 

i  =  1 
/ 

[ 

I  (f  v) 

3 

/ 

[H(1,  0),  L(p)J 

coercif 
op_principal 
[H(1,  0),  L(p) ] 

1 

*•**  hemicontlnu  (non  developpe) 

op_principsl 

[H(1,  0),  L(p)3 

•••»  borne  •••»  (non  developpe) 

°P_Princ*pal 
IH(1,  0),  L(p) 3 

monotone  (non  developpe) 

op_j>rincipal 
CH( 1 ,  0),  L(p)3 

aolution  •••* 

OKI.  0),  L(p)J 

»*>  fin 

TIm*  13181  Mac. 


Figure  6b.  Well-posedness  of  a  nonlinear  PDE. 
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Appendix:  Automatically  generated  FORTRAN  programs 

belman  .fortran 

subroutine  prodyn(nl,n2,epsimp,lmpmax,v,ro,u,eps,nmax) 
dimension  v(nl,n2),u(2,nl,n2) 

c  Resolution  de  1  equation  de  Bellman  dans  le  cas  ou: 
c  Les  parametres  sont 

c  L  etats-temps  est:  xl  x2 

c  La  dynamique  du  systeme  est  decrite  par  I  operateur 


C  Minu(  (0.5  -  Xl)**2  (1  -  x2)**2  X2**2  +  (1  -  Xl)  Xl  (1  -  X2)  x2 
c 

c  +  2  (l  -  x2)  x2  +  (l  -  Xl)**2  Xl**2  (0.5  -  x2)**2  +  2  (1  -  Xl)  Xl  +  u2**2  +  p2  U2 
c  +  Ul**2  +  pi  Ul  +  q2  +  ql  ) 
c 

c  ouv  designe  le  cout  optimal 

c  ou  pi  designe  sa  derivee  premiere  par  rapport  a  xi 

c  ou  qi  designe  sa  derivee  seconde  pan-  rapport  a  xi 

c  Le  probleme  est  statique 

c  Les  conditions  aux  limites  sont: 

c  x2  =  0  v=  xl  x2 

c  x2  =  1  v==  (l  -  xl)  (1  -  x2) 

c  xl  =  0  v==  xl  x2 

c  xl  =  1  v=  (1  -  xl)  (1  -  x2) 

c  Les  nombres  de  points  de  discretisation  sont:  nl  n2 
c  x2  =  1  correspond  a  12  =  n2 

c  x2  —  0  correspond  a  12  =  1 

c  xl  —  1  correspond  a  il  ==  nl 

c  xl  =  0  correspond  a  il  =  1 

c  Le  taux  d  actualisation  vaut:  1 

c  impmax  designe  le  nbre  maxi  d  iterations  du  systeme  implicite 
c  epsimp  designe  1  erreur  de  convergence  du  systeme  implicite 
c  ro  designe  le  pas  de  la  resolution  du  systeme  implicite 
c  par  une  methode  iterative 

c  Minimisation  par  la  methode  de  Newton  de  l’Hamiltonien 
c  L  inversion  de  la  Hessienne  est  faite  formellement 
c  nmax  designe  le  nombre  maxi  d  iteration  de  la  methode  de  Newton 
c  eps  designe  1  erreur  de  convergence  de  la  methode  de  Newton 

h2  ==  float(l)/(n2-l) 
hi  =  float(l)/(nl-l) 

U2  —  u(2,l,l) 

Ul  =  u(l,l,l) 
hih2  =  h2**2 
hihl  =  hl**2 
h22  =  2*h2 
h2l  =  2*hl 
nm2  =  n2-l 
nml  =  nl-1 
do  119  12  —  1  ,  n2  ,  l 
do  119  il  =  1  ,  nl  ,  l 
v(il,i2)  =  0.0 
119  continue 
imiter  =  1 
113  continue 
erimp  =  0 

do  ill  il  —  l  ,  nl  ,  l 
xl  ==  hi*(il-l) 
v(il,n2)  —  0 
v(il,l)  =  0 
111  continue 

do  109  12  =  2  ,  nm2  ,  1 

x2  -  -  Ii2*(i2-1) 
v(nl . i  —  l  —  0 
V(l,i2j  —  0 


110  continue 

do  109  il  =  2  ,  nml  ,  1 
xl  =  hl*(il-l) 

q2  —  (v(il,12+l)-2*v(il,i2)+v(il,i2-l))/hih2 
ql  =  (v(il  +  l,i2)-2*v(il,12)+v(il-l,i2))/hihl 
p2  =  (v(il,i2+l)-v(il,i2-l))/h22 
pi  =  (v(U+l,i2)-v(il-l,12))/h21 
niter  =  0 
wo  =  -1.0e+20 

101  continue 
niter  =  niter+i 

if  (  niter  -  nmax  )  102  ,  102  ,  103 

103  continue 
write(8,90l)il,i2 

901  format(’  newton  n  a  pas  converge’,  2  13) 
goto  104 

102  continue 

ui  =  -pl/2.0 
U2  =  -p2/2.0 

WW  =  (0.5-xl)**2*(l-x2)**2*x2**2+(l-xl)*xl*(l-x2)*x2+2*(l-x2)*x2+( 

1  l-xl)**2*xl**2*(0.5-x2)**2+2*(l-xl)*xl+u2**2+p2*u2+ul*»2+pl*ul  + 

2  q2+ql 

er  =  abs(ww-wo) 
if  (  er  -  eps  )  104  ,  104  ,  105 
105  continue 
wO  =  ww 
goto  101 

104  continue 
u(l,il,12)  =  ul 
u(2,il,i2)  =  u2 
w0  =  ww 

wo  =  wO-v(ll  ,12) 
vnew  =  ro*w0+v(il,l2) 
v(il,i2)  =  vnew 
erimp  =  abs(w0)+erimp 
109  continue 

imiter  =  imiter+1 

if  (  imiter  -  impmax  )  116  ,  115  ,  115 

116  continue 

if  (  epsimp  -  erimp  )  113  ,  112  ,  112 
115  continue 
write(8,907) 

907  format(’  schema  implicite  n  a  pas  converge’) 

112  continue 

do  117  il  =  1  ,  nl  ,  1 
do  117  i2  =  l  ,  n2  ,  l 
write(8,900)il,i2,v(il  ,12) 

900  format(’  v[’,  (13,’,’),  13,’]:’,  el4.7,'$’) 
wri  te(8, 902)11, i2,u(l,ll, 12), u(2,il,i2) 

902  format('  u[\  (13,’,'),  13, (el4,7,’,’),  el4.7,’]$’) 

117  continue 
return 
end 


c 

c  This  is  a  routine  that  solves  for  the  conditional  estimate  of 
c  x  at  t,  given  the  observations  y  up  to  time  t.  It  is  written 
c  by  Macsyma 
c 
c 

c  Program  Variables 
c 

c  xinput  real(4000) 

c  contains  the  actual  x  values  which  are  to  be  estimated  for 

c  times  at  instants  of  0.001  sec. 

c 

c  observ  real(4000) 

c  contains  the  sample  observations  to  be  filtered  again  at 

c  0.001  sec  time  intervals. 

c 

c  estim  real  (4000) 

c  contains  the  calculated  estimate  given  the  observations  up  to 

c  that  time  instant, 

c 

real  xinput(4002),  var(4002),  observ(4002),  estim(4002),  g(10),  ng 
1  (10),  xOxfx,  xlxfx,  x2xfx,  x3xrx,  x4xfx 

integer  i,  j,  k 
do  300  ii  —  1  ,  300  ,  1 

j=l 
k— 0 

size=0.00i 

xinput(0)=l 

observ(0)=l 

estim(0)=l 

do  100  i  =  1  ,  10  ,  1 
g(i)— o 
100  continue 
110  continue 

if(g(l)  .gt.  2.5)  goto  120 

call  nexty(xinput, observ,  1, size, J,k) 
deltay=observ(j)-observ(k) 
c 

sinhgl=sinh(l.O*g(l)) 

sinhg2=sinh(2.0*g(l)) 

coshgl=cosh(l.O*g(l)) 

coshg2=cosh(2.0*g(l)) 

ng(l)=size+g(l) 

c 

c 

ng(2)=coshg2*deltay+coshgl*deltay-1.0*deltay+g(2) 

c 

c 

ng(3)=-2.0*deltay*slnhg2+2.0*deltay*sinhgH-g(3) 

c 

c 

ng(4)=coshgl*deltay+g(4) 

c 

c 

ng(5)=g(5)-1.0*deltay*slnhgl 


c 

c 

ng(6)=coshg2*deltay-2.0*coshgl*deltay+deltay+g(6) 

c 

c 

cstl=-coshgl*sinhg2*size+coshg2*sinhgl*size+2.0*coshgl*sinhgl* 
1  size-sinhgl*size+coshgl*deltay**2*slnhg2 

cst2=-coshg2*deltay**2*sinhgl-2.0*coshgl*deltay**2*sinhgl+delt 
1  ay**2*sinhgl-2.0*g(2)*deltay*sinhgl-g(3)*coshgl*deltay 

ng(7)=cst2+cstl+g(7) 
c 
c 

cst3=-sinhgl*sinhg2*size4-sinhgl**2*size-coshgl*coshg2*size+2.0 
1  *coshgl**2*size-coshgl*size 

cst4=deltay**2*sinhgl*sinhg2-deltay**2*slnhgl**2-g(3)*deltay*s 
l  inhgl+coshgl*coshg2*deltay  **2-2.0*coshgl  **2*deltay**2 

cst5=coshgl*deltay**2+2.0*g(4)*coshg2*deltay-4.0*g(4)*coshgl*d 
1  eltay+2.0*g(4)*deltay+g(8) 

ng(8)=cst5+cst4+cst3 
c 
c 

ng(9)=0.5*coshgi*sinhgl*slze-0.5*coshgl*deltay**2*sinhgl-1.0*g 
1  (4)*deltay*sinhgl+g(9) 

c 

c 

cst6=-g(4)*sinhgl*sinhg2*slze+g(5)*coshgl*sinhg2*size+g(4)*sin 
X  hgl**2*size+g(2)*sinhgl**2*size-g(5)*coshg2*sinhgl*size 

cst7=-2.0*g(5)*coshgl*sinhgl*size+g(3)*coshgl*sinhgl*size+g(5) 

1  *sinhgl*size-g(4)*coshgl*coshg2*size+2.0*g(4)*coshgl**2*siz 

2  e 

cst8=-g(4)*coshgl*size+g(4)*deltay**2*slnhgl*sinhg2-g(5)*coshg 
1  I*deltay**2*sinhg2-deltay*sinhg2-g(4)*deltay**2*sinhgl**2 

cst9=-g(2)*deltay**2*sinhgl**2+g(5)*coshg2*deltay**2*sinhgl+2. 

1  0*g(5)*coshgl*deltay**2*sinhgl-g(3)*coshgl*deltay**2*sinhgl 

2  -g(5)*deltay**2*sinhgl 

cstl0=2.0*g(2)*g(5)*deltay*sinhgl-g(3)*g(4)*deltay*sinhgl+delt 

1  ay*sinhgl+g(4)*coshgl*coshg2*deltay**2-2.0*g(4)*coshgl**2*d 

2  eltay**2 

cstll=g(4)*coshgl*deltay**2+g(4)**2*coshg2*deltay+g(3)*g(5)*co 
1  shgl*deltay-2.0*g(4)**2*coshgl*deltay+g(4)**2*deltay 

ng(10)=cst9+cst8+cst74-cst6+cstll+cstl0+g(10) 
c 
c 

do  105  i  =  1  ,  10 , 1 
g(i)=ng(i) 

105  continue 
c 

xOxfx=xOfx(g) 

xlxfx=xlfx(g) 

x2xrx=x2fx(g) 

x3xfx=x3fx(g) 

x4xfx=x4fx(g) 

sinhgl=slnh(l.0*g(l)) 

slnhg2=slnh(2.0*g(l)) 

sinhg3=sinh(3 ,0*g(  1 )) 


coshgl=cosh(1.0*g(l)) 

coshg3=cosh(3.0*g(l)) 

coshg2=cosh(2.0*g(l)) 

cstl2=-0.25*coshgl**2*sinhg3*x3xfx/sinhg  1**2-0.083333333333333 

1  33*sinhg3*x3xfx-0.75*sinhgl*x3xfx+0.25*coshgl*coshg3*x3xfx/ 

2  sinhgl+1.5*coshgl**2*x3xfx/sinhgl 
cstl3=-coshgl*x3xfx/sinhgl-|-0.08333333333333333*coshgl**3*coshg 

1  3*x3xfx/sinhgl**3-0.75*coshgl**4*x3xfx/sinhgl**3-t-0.66666666 

2  66666667  *coshgl**3*x3xfx/sinhgl*  *3-0.5*  g(5)*coshgl*sinhg3*X 

3  2xfx/sinhgl**2 

csti4=0.5*coshgl*sinhg3*x2xfx/sinhgl**2-g(6)*coshgl*sinhg2*x2x 

1  fx/sinhgi-g(2)*coshgl*sinhg2*x2xfx/sinhgl+0.5*g(3)*coshgl** 

2  2*sinhg2*x2xfx/sinhgl**2+0.5*g(3)*sinhg2*x2xfx 
cstl5=0.25*g(5)*coshg3*x2xfx/sinhgl-0.25*coshg3*x2xfx/sinhgl-g 

1  (3)*coshgl*coshg2*x2xfx/sinhgl+2.25*g(5)*coshgl*x2xfx/sinhg 

2  l-2.25*coshgl*x2xfx/sinhgl 

cstl6— -g(5)*x2xfx/sinhgl+x2xfx/sinhgl+0.5*g(6)*coshgl**2*coshg 

1  2*x2xfx/sinhgl**2+0.5*g(2)*coshgl**2*coshg2*x2xfx/sinhgl**2 

2  +0.5*g(6)*coshgl**2*x2xfx/sinhgl**2 
cstl7=-0.5*g(2)*coshgl**2*x2xfx/slnhgl**2+0.25*g(5)*coshgl**2* 

1  coshg3*x2xfx/sinhgl**3-0.25*coshgl**2*coshg3*x2xfx/sinhgl** 

2  3-2.25*g(5)*coshgl**3*x2xfx/sinhgl**3+2.25*coshgl**3*x2xfx/ 

3  sinhgl**3 

cstl8=2.0*g(5)*coshgl**2*x2xfx/sinhgl**3-2.0*coshgl**2*x2xfx/s 

1  inhgl**3+0.5*g(6)*coshg2*x2xfx+0.5*g(2)*coshg2*x2xfx-0.5*g( 

2  6)*x2xfx 

csti9=0.5*g(2)*x2xfx-t-0.5*coshgl*sinhg3*xlxfx/slnhgl-0.25*g(5)* 

1  *2*sinhg3*xlxfx/sinhgl**2+0.5*g(5)*sinhg3*xlxftc/sinhgl**2-0 

2  .25*sinhg3*xlxfx/sinhgl**2 
cst20=-g(5)*g(6)*sinhg2*xlxfx/sinhgl+g(6)*sinhg2*xlxfx/sinhgl- 

1  g(2)*g(5)*sinhg2*xlxfx/sinhgl+g(2)*sinhg2*xlxfx/sinhgl+g(3) 

2  *g(5)*coshgl*sinhg2*xlxfx/slnhgl**2 

cst2l— -g(3)*coshgl*sinhg2*xlxfx/sinhgl**2+g(8)*sinhgl*xlxfx-2. 

1  0*g(4)*g(6)*sinhgl*xlxfx-g(3)*g(5)*coshg2*xlxfx/sinhgl+g(3) 

2  *coshg2*xlxfx/sinhgl 

cst22=-g(8)*coshgl**2*xlxfx/sinhgl+2.0*g(4)*g(6)*coshgl**2*xlx 

1  fx/sinhgl+0.75*g(5)**2*xlxfx/sinhgl-1.5*g(5)*xlxfx/sinhgl+0 

2  .75*xlxfx/sinhgl 

cst23=-0.25*coshgl**2*coshg3*xlxfx/sinhgl**2+g(5)*g(6)*coshgl* 

1  coshg2*xlxfx/sinhgl**2-g(6)*coshgl*coshg2*xlxfx/sinhgl**2+g 

2  (2)*g(5)*coshgl*coshg2*xlxfx/sinhgl**2-g(2)*coshgl*coshg2*X 

3  lxfx/sinhgl**2 

cst24— 2.25*coshgl**3*xlxfx/sinhgl**2-2.0*coshgl**2*xlxfx/sinhg 

1  l**2+g(5)*g(6)*coshgl*xlxfx/sinhgl**2-g(6)*coshgl*xlxfx/sin 

2  hgl**2-g(2)*g(5)*coshgl*xlxfx/slnhgl**2 
cst25=g(2)*coshgl*xlxfx/sinhgl**2+0.25*g(5)**2*coshgl*coshg3*x 

1  lxfx/sinhgl**3-0.6*g(5)*coshgl*coshg3*xlxfx/sinhgl**3+0.25* 

2  coshgl*coshg3*xlxfx/sinhgl**3-2.25*g(5)**2*coshgl**2*xlxfx/ 

3  sinhgl**3 

cst26==4.5*g(5)*coshgl**2*xlxfx/sinhgl**3-2.25*coshgl**2*xlxfx/ 

1  slnhgl**3+2.0*g(5)**2*coshgl*xlxfx/slnhgl**3-4.0*g(5)*coshg 

2  l*xlxfx/sinhgl**34-2.0*coshgl*xlxfx/8inhgl**3 
cst27=-0.25*coshg3*xixfx-2.25*coshgl*xlxrx-i-xlxrx-l-0.26*g(6)*sin 

1  hg3*x0xfx/slnhgl-0.25*slnhg3*x0xfx/slnhgl 

cst28=-0.5*g(3)*coshgl*sinhg2*x0xfx/slnhgl-t-0.5*g(3)*g(5)**2*si 


1  nhg2*x0xfx/sinhgl**2-g(3)*g(5)*sinhg2*x0xfx/sinhgl**2+0.5*g 

2  (3)*sinhg2*x0xfx/sinhgl**2+0.5*g(6)*sinhg2*x0xfx 
cst29=0.5*g(2)*sinhg2*x0xfx-0.5*g(6)*coshgl*coshg2*x0xfx/sinhg 

1  l-0.5*g(2)*coshgl  *coshg2*xOxfx/sinhgl-g(5)*g(8)*coshgl*xOxf 

2  x/sinhgl+g(8)*coshgl*xOxfx/sinhgl 
cst30=2.0*g(4)*g(5)*g(6)*coshgl*x0xfx/sinhgl-2.0*g(4)*g(6)*cos 

1  hgl*x0xfx/sinhgl-0.5*g(6)*coshgl*x0xfx/sinhgl+0.5*g(2)*cosh 

2  gl*x0xfx/sinhgl-0.25*g(5)*coshgl*coshg3*x0xfx/sinhgl**2 
cst31=0.25*coshgl*coshg3*x0xfx/sinhgl**2+0.5*g(5)**2*g(6)*cosh 

1  g2*x0xfx/sinhgl**2-g(5)*g(6)*coshg2*x0xfx/sinhgl**2+0.5*g(6 

2  )*coshg2*x0xfx/sinhgl**2+0.5*g(2)*g(5)**2*coshg2*x0xfx/sinh 

3  gl**2 

cst32=-g(2)*g(5)*coshg2*x0xfx/slnhgl**2+0.5*g(2)*coshg2*x0xfx/ 

1  sinhgl**2+2.25*g(5)*coshgl**2*xOxfx/sinhgl**2-2.25*coshgl** 

2  2*x0xfx/sinhgl**2-2.0*g(5)*coshgl*x0xfx/sinhgl**2 
cst33=2.0*coshgl*x0xfx/slnhgl**2+0.5*g(5)**2*g(6)*x0xfx/sinhgl 

X  **2-g(5)*g(6)*x0xfx/sinhgl**2+0.5*g(6)*x0xfx/sinhgl**2-0.5* 

2  g(2)*g(5)**2*xOxfx/sinhgl**2 

cst34=g(2)*g(5)*x0xfx/sinhgl**2-0.5*g(2)*x0xfx/sinhgl  **2+0.083 

1  33333333333333  *g(5)**3*coshg3*x0xfx/slnhgl**3-0.25*g(5)*  *2* 

2  coshg3*x0xfx/sinhgl**3+0.25*g(5)*coshg3*x0xfx/sinhgl**3 
cst35=-0.08333333333333333*coshg3*x0xfx/sinhgl**3-0.75*g(5)**3 

1  *coshgl*x0xfx/slnhgl**3+2.25*g(5)**2*coshgl*x0xfx/sinhgl**3 

2  -2.25*g(5)*coshgl*x0xfx/sinhgl**3+0.75*coshgl*x0xfx/sinhgl* 

3  *3 

cst36=0. 6666666666666667  *g(5)**3*x0xfx/slnhgl  **3-2. 0*g(5)**2*X 

1  Oxfx/sinhgl*  *3+ 2. 0*g(5)*x0xfx/sinhgl*  *3-0. 6666666666666667* 

2  x0xfx/slnhgl**3+0.5*g(3)*coshg2*x0xfx 
Cst37=g(l0)*x0xfx-g(4)*g(8)*x0xfx+g(7)*x0xfx+g(4)**2*g(6)*x0xf 

1  x-0.75*g(5)*x0xrx 

Ul=-0.5*g(3)*x0xfx+0.75*x0xfx+cst37+cst36+cst35+cst34+cst33+cs 

1  t32+cst31+cst30+cst29+cst28+cst27+cst26+cst25+cst24+cst23+c 

2  st22+cst21+cst20+cstl9+cstl8+cstl7+cstl6+cstl5+cstl4+cstl3+ 

3  cstl2 

cst38=-0.25*coshgl**2*slnhg3*x4xrx/sinhgl**2-0.083333333333333 

1  33*slnhg3*x4xfx-0.75*sinhgl*x4xfx+0.25*coshgl*coshg3*x4xfx/ 

2  sinhgl+1.5*coshgl**2*x4xfx/slnhgl 
cst39=-coshgl*x4xfx/slnhgl+0.08333333333333333*coshgl**3*coshg 

1  3*x4xfx/sinhgl*  *3-0 .75*coshgl**4*x4xfx/slnhgl*  *3+0. 66666666 

2  66666667*coshgl**3*x4xfx/sinhgl**3-0.5*g(5)*coshgl*sinhg3*x 

3  3xfx/sinhgl**2 

cst40=0.5*coshgl*sinhg3*x3xfx/sinhgl**2-g(6)*coshgl*sinhg2*x3x 
X  fx/sinhgl-g(2)*coshgl*sinhg2*x3xfx/sinhgl+0.5*g(3)*coshgl** 

2  2*sinhg2*x3xfx/sinhgl**2+0.5*g(3)*sinhg2*x3xfx 

cst41=0.25*g(5)*cosXig3*x3xrx/slnhgl-0.25*coshg3*x3xfx/sinhgl-g 

1  (3)*coshgl*coshg2*x3xfx/sinhgl+2.25*g(5)*coshgl*x3xfx/sinhg 

2  l-2.25*coshgl*x3xfx/sinhgl 
cst42=-g(5)*x3xfx/sinhgl+x3xfx/sinhgl+0.5*g(6)*coshgl**2*coshg 

1  2*x3xfx/slnhgl**2+0.5*g(2)*coshgl**2*coshg2*x3xfx/8inhgl**2 

2  +0.6*g(6)*coshgl**2*x3xfx/slnhgl**2 

cst43= -0.5*g(2)*coshgl**2*x3xfx/sinhgl**2+0.25*g(6)*COShgX**2* 

1  cO'lig3*x3Xfx/sinhgl**3-0.25*coshgl**2*coshg3*x3xfx/sinhgl** 

2  3-2  25*g(5)*coshgl**3*x3xfx/sinhgX**3+2.25*co8hgl**3*x3xfx/ 

3  feinligl**3 


cst44=2.0*g(5)*coshgl**2*x3xfx/sinhgl**3-2.0*coshgl**2*x3xfx/s 

1  inhgl**3+0.5*g(6)*coshg2*x3xfx-f0.5*g(2)*coshg2*x3xfx-0.5*g( 

2  6)*x3xfx 

cst45=0.5*g(2)*x3xfx+0.5*coshgl*sinhg3*x2xfx/sinhgl-0.25*g(5)* 

1  *2*sinhg3*x2xfx/sinhgl**2+0.5*g(5)*sinhg3*x2xfx/sinhgl**2-0 

2  ,25*sinhg3*x2xfx/sinhgl**2 
cst46=-g(5)*g(6)*sinhg2*x2xfx/sinhgl+g(6)*sinhg2*x2xfx/sinhgl- 

1  g(2)*g(5)*sinhg2*x2xfx/slnhgl+g(2)*sinhg2*x2xfx/sinhgl+g(3) 

2  *g(5)*coshgl*sinhg2*x2xfx/sinhgl**2 
cst47=-g(3)*coshgl*sinhg2*x2xfx/sinhgl**2+g(8)*sinhgl*x2xfx-2. 

X  0*g(4)*g(6)*sinhgl*x2xfx-g(3)*g(5)*coshg2*x2xfx/sinhgH-g(3) 

2  *coshg2*x2xfx/sinhgl 

cst48=-g(8)*coshgl**2*x2xfx/sinhgH-2.0*g(4)*g(6)*coshgl**2*x2x 
X  fx/sinhgl+0.75*g(5)**2*x2xfx/sinhgl-1.5*g(5)*x2xfx/sinhgl+0 

2  .75*x2xfx/sinhgl 

cst49=-0.25*coshgl**2*coshg3*x2xfx/sinhgl**2+g(5)*g(6)*coshgl* 

1  coshg2*x2xfx/sinhgl**2-g(6)*coshgl*coshg2*x2xfX/sinhgl**2+g 

2  (2)*g(5)*coshgl*coshg2*x2xfx/sinhgl**2-g(2)*coshgl*coshg2*x 

3  2xfx/sinhgl**2 

cst50=2.25*coshgl**3*x2xfx/sinhgl**2-2.0*coshgl**2*x2xfx/sinhg 

1  l**2+g(5)*g(6)*coshgl*x2xfx/sinhgl**2-g(6)*coshgl*x2xfx/sin 

2  hgl**2-g(2)*g(5)*coshgi*x2xfx/sinhgi**2 
cst51=g(2)*coshgl*x2xfx/slnhgl**2+0.25*g(5)**2*coshgl*coshg3*x 

1  2xfx/sinhgl**3-0.5*g(5)*coshgl*coshg3*x2xfx/sinhgl**3+0.25* 

2  coshgl*coshg3*x2xfx/sinhgl**3-2.25*g(5)**2*coshgl**2*x2xfx/ 

3  sinhgl**3 

cst52=4.5*g(5)*coshgl**2*x2xfx/sinhgl**3-2.25*coshgl**2*x2xfx/ 

1  sinhgl**3+2.0*g(5)**2*coshgl*x2xfx/sinhgl**3-4.0*g(5)*coshg 

2  I*x2xfx/sinhgl**3+2.0*coshgl*x2xfx/sinhgl**3 
cst53=-0.25*coshg3*x2xfx-2.25*coshgl*x2xfx+x2xfx+0.25*g(5)*sin 

X  hg3*xlxfx/sinhgl-0.25*sinhg3*xlxfx/sinhgl 

cst54=-0.5*g(3)*coshgl*sinhg2*xlxfx/sinhgl+0.5*g(3)*g(5)**2*si 

1  nhg2*xlxfx/slnXigl**2-g(3)*g(5)*sinhg2*xlxfX/sinhgl**2+0.5*g 

2  (3)*sinhg2*xlxfx/sinhgl**2+0.5*g(6)*sinhg2*xlxfx 
cst55=0.5*g(2)*sinhg2*xlxrx-0.5*g(6)*coshgl*coshg2*xlxfX/sinhg 

1  l-0.5*g(2)*coshgl*coshg2*xlxfx/sinhgl-g(5)*g(8)*coshgl*xlxf 

2  x/sinhgl+g(8)*coshgl*xlxfx/sinhgl 
cst56=2.0*g(4)*g(5)*g(6)*coshgl*xlxfx/sinhgl-2.0*g(4)*g(6)*cos 

1  hgl*xlxfx/sinhgl-0.5*g(6)*coshgl*xlxfx/sinhgl+0.5*g(2)*cosh 

2  gl*xlxfx/sinhgl-0.25*g(5)*coshgl*coshg3*xlxfx/sinhgl**2 
cst57=0.25*coshgl*coshg3*xlxrx/sinhgl**2+0.5*g(5)**2*g(6)*cosh 

1  g2*xlxfx/sinhgl**2-g(5)*g(6)*coshg2*xlxfx/sinhgl**2+0.5*g(6 

2  )*coshg2*xlxfx/slnhgl**2+0.5*g(2)*g(5)**2*coshg2*xlxfx/sinh 

3  gl**2 

cst58=-g(2)*g(5)*coshg2*xlxfx/sinhgl**2+0.5*g(2)*coshg2*xlxfx/ 

1  sinhgl**2+2.25*g(5)*coshgl**2*xlxfx/sinhgl**2-2.25*coshgl** 

2  2*xlxfx/sinhgl**2-2.0*g(5)*coshgl*xlxfx/sinhgl**2 
cst59==2.0*coshgl*xlxfx/sinhgl**2+0.5*g(5)**2*g(6)*xlxfx/sinhgl 

1  **2-g(5)*g(6)*xlxfx/sinhgl**2+0.5*g(6)*xlxfx/sinhgl**2-0.5* 

2  g(2)*g(5)**2*xlxfx/slnhgl**2 
cst60=g(2’l*g(5)*xlxfx/sinhgl  **2-0. 5*g(2)*xlxfx/sinhgl  **2+0.083 

1  3333,' ’,:!.'..i333333*g(5)**3*coshg3*xlxrx/6inhgl  **3-0. 25*g(5)**2* 

2  coslig:!'  \  lxfx/slnhgl**3+0.25*g(5)*coshg3*xlxrx/slnhgl**3 
cst61  =  -0  O8333333333333333*coshg3*xlxfx/sinhgl**3-0.75*g(6)**3 

1  *coshgi  *\\lxfx/sinhgl**3+2.25*g{5)**2*coshgl*xlxfx/sinhgl**3 


2 

3 


-2.25*g(5)*coshgl*xlxfx/sinhgl**3+0.76*coshgl*xlxfx/sinhgl* 

*3 

cst62— 0.6666666666666667*  g(5)**3*xlxfx/sinhgl  **3-2. 0*g(5)**2*X 

1  lxfx/sinhgl**3+2.0*g(5)*xlxfx/sinhgl**3-0. 6666666666666667* 

2  xlxfx/sinhgl**34-0.5*g(3)*coshg2*xlxfx 
cst63=g(l0)*xlxfx-g(4)*g(8)*xlxfx+g(7)*xlxfx+g(4)**2*g(6)*xlxf 

1  x-0.75*g(5)*xlxfx 

ulx=-0.5*g(3)*xlxfx-)-0.75*xlxfx4-cst63+cst62-l-cst61-|-cst60-l-cst59-|-c 

1  st58-fcst57+cst56-|-cst55+cst54-i-cst53-|-cst52-)-cst51-fcst50+ cst49+ 

2  CSt48-)-CSt47-fCSt46-|-CSt45+CSt44+CSt43+cst42+CSt4H-CSt40-)-cst39 

3  +CSt38 
C 

estim(  j  )=(ulx-u  l*xlxfx/x0xfx)/x0xfx+x  lxfx/xOxfX 

var(j)=(xinput(j)-estim(j))**2-t-var(j) 

k=j 

J=j+l 

go  to  1X0 
120  continue 
300  continue 

do  400  i  =  1  ,  4000 , 1 
var(i)=var(i)/300.0 
y=i*size 
print*,y,var(i) 

400  continue 
stop 
end 

real  function  xOfx(g) 
real  g(l) 

xOfx=l  /sqrt(cosh(g(l))) 

return 

end 

real  function  xlfx(g) 
real  g(x) 

xlfx=(l-g(5))*cosh(g(l))**((-3.0)/2.0) 

return 

end 

real  function  x2fx(g) 
real  g(l) 

x2fx==cosh(g(l))**((-5.0)/2.0)*(cosh(g(l))*sinh(g(l))+g(5)**2-2*g(5 

1  )+l) 
return 
end 

real  function  x3fx(g) 
real  g(l) 

x3fx=(l-g(5))*cosh(g(l))**((-7.0)/2.0)*(3*cosh(g(l))*sinh(g(l))+g( 

1  5)**2-2*g(6)+l) 

return 
end 

real  function  x4fx(g) 
real  g{l) 


x4fx=cosh(g(l))**((-0.O)/2.O)*(3*cosh(g(l))**2+sinh(g(l))**2+6*g(5 

1  )**2*cosh(g(l))*sinh(g(l))-12*g(5)*cosh(g(l))*sinh(g(l))+6*cosh 

2  (g(l))*slnh{g(l))+g(5)**4-4*g(5)**3+6*g(5)**2-4*g(5)+l) 
return 

end 

subroutine  nexty(xinput,observ,e,8lze,j,k) 
real  xlnput(l),  observ(l),  e,  size 
integer  j,  k 

xinput(j)=grand(0)*sqrt(slze)+xlnput(k) 

observ(j)=(e*xinput(j)**2+xinput(j))*slze+grand(0)*sqrt(slze)+obse 

1  rv(k) 
return 
end 


